#clean up
rm(list=ls())

#set seed
set.seed(715130425)

#load libraries
library(akima)
library(geoR)
library(scatterplot3d)  
library(ggmap)
library(lattice)
library(maps)
library(spdep)
library(doParallel)
library(foreach)

#options
options(scipen=8)

#Register Number of Processors for Parallel Computing
cores<-detectCores()-1; cores
#cores<-7; cores
registerDoParallel(cores)
getDoParWorkers()


###BAYESIAN MODEL USING BOOTSTRAPPED SAMPLES###
#number of iterations
K<-100
getDoParWorkers()

###start parallel loop###
out<-foreach(k=1:K, .combine=list, .inorder=FALSE, .errorhandling='remove') %dopar%{
#out<-foreach(k=1:K, .combine=list, .inorder=TRUE, .errorhandling='stop') %dopar%{
	
	#load subset of original data
	subset.data<-read.csv(paste('bootstrap/train',k,'.csv',sep=""))

	#load 5 percent of kriged points
	subset.covariates<-read.csv(paste('bootstrap/test',k,'.csv',sep=""))

	#number of simulated observations
	num.sims<-dim(subset.covariates)[1]

	#estimate model on subset
	ideol.bayes.local <- krige.bayes(coords=cbind(subset.data$eastings,subset.data$northings), data=subset.data$ideology, locations=NULL,
		model = model.control(trend.d =~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
		subset.data$inc14+subset.data$cathOrth+subset.data$consRelig+subset.data$musJew+
		subset.data$mainline+subset.data$rural+subset.data$ownership+as.factor(subset.data$empstat)+
		subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2),
		cov.model = "gaussian"),
		prior = prior.control(beta.prior="flat", tausq.rel.prior="uniform", tausq.rel.discrete=seq(from=6,to=8,by=.05)))

	#iterative output of posterior
	postSample<-ideol.bayes.local$posterior$sample
	write.csv(postSample,paste('bootstrap/sample',k,'.csv',sep=""),row.names=F)

	#simple kriging by iteration
	#create output matrices
	predictions<-matrix(NA,nrow=1000,ncol=num.sims)
	variances<-matrix(NA,nrow=1000,ncol=num.sims)
	
	#loop through each iteration
	for(i in 1:1000){
		kriges<-krige.conv(coords=cbind(subset.data$eastings,subset.data$northings), data=subset.data$ideology, locations=cbind(subset.covariates$eastings,subset.covariates$northings),
			krige=krige.control(type.krige='sk', trend.d=~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
		subset.data$inc14+subset.data$cathOrth+subset.data$consRelig+subset.data$musJew+
		subset.data$mainline+subset.data$rural+subset.data$ownership+as.factor(subset.data$empstat)+
		subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2),
			trend.l=~subset.covariates$age*subset.covariates$educ+as.factor(subset.covariates$race)*subset.covariates$female+
		subset.covariates$inc14+subset.covariates$cathOrth+subset.covariates$consRelig+subset.covariates$musJew+
		subset.covariates$main+subset.covariates$rural+subset.covariates$ownership+as.factor(subset.covariates$empstat)+
		subset.covariates$eastings*subset.covariates$northings+I(subset.covariates$eastings^2)+I(subset.covariates$northings^2),
			beta=as.numeric(postSample[i,1:23]), cov.pars=c(postSample$sigmasq[i],postSample$phi[i]), 
			nugget=I(postSample$sigmasq[i]*postSample$tausq.rel[i]), cov.model='gaussian'))
		predictions[i,]<-kriges$predict
		variances[i,]<-kriges$krige.var
	}
	
	#mean for each location
	predictive.mean<-apply(predictions,2,mean)
	
	#variance for each location
	between<-apply(predictions,2,var)
	within<-apply(variances,2,mean)
	predictive.variance<-within+((num.sims+1)/num.sims)*between
	
	#create matrix of kriged points
	forecasts<-cbind(subset.covariates,predictive.mean,predictive.variance)
	#forecasts<-cbind(subset.covariates[,c(1,9:14,24)],predictive.mean,predictive.variance)

	#iterative output of forecasts
	write.csv(forecasts,paste('bootstrap/forecast',k,'.csv',sep=""),row.names=F)

rm(subset.data,subset.covariates,postSample,forecasts,num.sims)
}



